/*  sdfexpds.cpp  by K.Tsuru  */
// function ID 3313  DRADIX
/**************************************************************************
SDouble class
It provides the exponential function exp(x).
From the condition "exp(x) < DRADIX^DRADIX_EXP_MAX" the range of 'x' must be
within
 |x| < DFIGURES*log(10)*DRADIX_EXP_MAX = 301759.17...

[Algorithm of ExpDS()](divide and series method)
step 1 :It divides 'x' into an integral part 'n' and a decimal one 'd'.
    x = n + d( n is integer and |d|<=0.5).
step 2 :It reduces 'd' by taking M=2^m and f = d/M.
step 3 :It evaluates y = e^f by the series.
    In this evaluation it divides 'f' such as
        f=S(short decimal)+L(long but very small decimal)
    and calculates y = (e^S)*(e^L).This processing uses the fact that
    the series(products) of short decimal can be rapidly evaluated.
    Using this it seems that step 2 is not necessary,but the calculation
    of e^S costs much time.
------------------------------------------------------------
The method in which the reduction method is used only to e^S and calculates
(e^(S/M))^M *(e^L) costs the same or more time.
------------------------------------------------------------
step 4 :It restores the reduction by e^d=y^M.
step 5 :It returns the value e^x = (e^n)*(e^d).
***************************************************************************/

#ifndef SN_H
#include "sn.h"
#endif

static const char* const func = "ExpDS";
//It uses the series when the number of figures is less than this value.
static const uint seriesFig = 15u;

/*********************************************
sub-function
From 'dp' it determines 'rPow2' and 'upPrec'.
**********************************************/
static void GetReduceMethod(const SDouble& dp, int* rPow2, uint* upPrec){
    uint fig = dp.Last() - dp.First() +1u;  //figure number of decimal part
    int k, exp = dp.NetRdxExp();
//When even if 'fig' is small but the effective figures is large,it takes large 'k'.
    if(dp.MaxSize() >= 512u) k = 32;
    else if(fig < 128u) k = 8;
    else if(fig < 256u) k = 16;
    else k = 32;
    ulong f = (ulong)dp(dp.First());    // f<DRADIX<2^14 = 16384
//Notice : ulong=32 bits, the value "f>>32" cannot be evaluated.
    if(!exp && (k < 32) && (f >> k) ){  //It makes the reduced value less than 1/DRADIX.
    //When k==8 only, it comes here because f/(2^16)=0.
        k *= 2;
    }
//It gives the number of figures which temporally raising up precision.
//That is the number of figures of 2^rPow2. The cancelation occures by this figures.
    *upPrec = uint( k*log10(2.0)/DFIGURES ) + 1u;
    *rPow2 = k;
}
/****************************
sub-function
It sets "result=exp(decPart)".
*****************************/
static void ExpLongDecPart(SDouble& decPart, SDouble& result){
    RealSize C;
    int j, rPow2;
    uint upPrec;
    GetReduceMethod(decPart, &rPow2, &upPrec);
    uint up = result.ProperUpPrec(upPrec);
    //It temporally raises up the precision if possible.
    if(up) C.SetEffFig(C.EffFigures() + up, C.TEMP_EXTEND);
//step 2: It devides "decPart" by "2^rPow2".
    ulong dv;
    j = rPow2;
    if(j <= 16){
        dv = 1uL << j;  // 2^16 = 65536 < SlOpMaxValue()
        decPart = DsDiv(decPart, dv);
    } else {
        dv = 1uL << 16;     j /= 16;
    // 2^rPow2 = (2^16)^(rPow2/16) :It devides 'j' times by 2^16 .
        while(j > 0){   decPart = DsDiv(decPart, dv);   j--;    }
    }
/*
It uses the fact that the series of short decimal can be rapidly
evaluated.
 e^d = e^s*e^(d-s) : 's' is a short decimal and |d-s|<<1.0.
When "EffFigures()/4=500",it becomes about twice faster.
*/
    SDouble rr, sd, one(1.0);
    uint tk = seriesFig, df;
    df = decPart.Last() - decPart.First() +1u;
//If the remaining figures is small it unifies.
    if(df/tk < 2u) tk = df;
// step 3.1 :It takes out upper 'tk& figures and splits 'd' into 's' and 'd-s'.
    sd = decPart.TakeOutFigures(tk);
    decPart -= sd;
// step 3.2:evaluation of exp(s)
    rr = ExpSeries(sd);     // sd : short decimal
#ifndef NDEBUG
if(decPart.Sign()) assert(decPart.NetRdxExp() <= -(int)seriesFig);
#endif
// step 3.3 :evaluation of exp(d-s)
    if(decPart.Sign()){
        result = ExpSeries(decPart);    // decPart :long but very small decimal fraction
        result *= rr;
    } else result = rr;
// step 4 :It restores the reduction. result = result^p, p = 2^rPow2
    j = rPow2;
    while(j){   //If the value of 'rPow2' is taken too large it costs much time here.
        result *= result;   j--;    //(...(result^2)^2....^2)^2
    }
    if(up){
        C.SetEffFig(0);     result.Reform(3301);    //It shortens within the effective figures.
    }
}

/********************
main body of ExpDS(x)
*********************/
static const double xMax = (double)DFIGURES*M_LN10*(double)DRADIX_EXP_MAX;
     // = 301759.17...
SDouble ExpDS(const SDouble& x){
    SDouble one(1.0);
    if( x.IsMLT(one) ) return (one + x);    // x = 0 or |x| << 1.0, x*x = 0;

    SDouble decPart, intPart;
// step 1 :It splits into an integral part and decimal one.
    decPart = Modf(x, intPart);
    double iX = doubleD(intPart, 0), dX = doubleD(decPart, 0);
// check the condition "exp(x) < DRADIX^DRADIX_EXP_MAX".
#ifndef NDEBUG
assert(xMax < (double)LONG_MAX);
#endif
    if( fabs(iX) > xMax){   // exp(x) > D^(DRADIX_EXP_MAX)
        if(iX > 0.0) x.SetError(x.OVERFLOW_ERR, func, 3301);
        x.SetError(x.UNDERFLOW_ERR, func, -3301);
        return 0.0;
    }
    long ePow = iX; // iX < xMax < LONG_MAX
//It reduces as |decPart| <= 0.5.
// example : x = -1.8 --> intPart = -2, decPart = 0.2
    if(dX > 0.5){
        ePow++;     decPart -= one;
    } else if(dX < -0.5){
        ePow--;     decPart += one;
    }

    SDouble ePowInt;
    if(ePow) {
        ePowInt = Expower(ePow);    //exp(ePow)
        if(decPart.Sign(3301) == 0) return ePowInt; // x is an integer.
    }
//evaluation of exp(decPart)
    SDouble result;
    uint df = decPart.Last() - decPart.First() +1u;//number of figures of decimal part
// seriesFig :If the number of figures is less than this value it uses the raw series.
    if( (df >= seriesFig) || (decPart.NetRdxExp() >= 0) ){  //long decimal
        ExpLongDecPart(decPart, result);
    } else {//If 'x' is short decimal and less than 1/DRADIX it is faster not to use the reduction.
            //It is better not to use the reduction method of x/(2^p) too.
        result = ExpSeries(decPart);
    }
// step 5 :multply by exp(n)
    if(ePow) result *= ePowInt;     // includes Reform()
    return result;
}